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ABSTRACT 

The  feasibility  of  applying  the  principles  of  matched 
field  processing  to  ocean  acoustic  tomography  were  studied 
under  various  conditions  of  ambient  noise.  Several  likeli- 
hood estimators  were  examined  (e.g.,  Bucker,  Bartlett, 
Maximum  Likelihood,  etc.).  Simulations  were  initially 
conducted  for  the  simple  case  wherein  only  one  parameter  of 
the  medium  was  unknown  (e.g.,  SOFAR  axis  depth,  surface 
sound  speed,  position  of  a  single  acoustic  front) .  The 
method  was  then  applied  to  the  more  realistic  problem  of 
locating  the  boundaries  of  an  eddy  in  the  ocean.  For 
moderate  signal-to-noise  ratios,  all  the  estimators  were 
shown  to  be  able  to  solve  the  problem,  albeit  with  different 
efficiencies.  For  low  signal-to-noise  ratios,  the  MLM 
scheme  proved  to  be  the  most  reliable  especially  when  a 
highly  correlated  ambient  noise  was  present.  In  all  cases, 
computer  simulations;  illustrated  that  mismatching  may  occur 
when  the  parameterization  of  the  medium  is  poorly 
approximated.  Mismatching  leads  to  a  decrease  in  the 
efficiency  of  the  estimators  but  it  may  be  still  possible  to 
correctly  estimate  the  environmental  characteristics. 
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I.   INTRODUCTION 

As  sound  waves  propagate  through  the  ocean,  the  complex 
acoustic  pressure  field  which  is  generated  by  the  source 
depends  mainly  on  the  path  followed  by  the  acoustic  rays  and 
the  sound  speed  along  this  particular  path.  Due  to  this 
close  relationship  between  the  sound  speed  field  and  the 
acoustic  pressure  field,  an  attempt  may  be  made  to  estimate 
the  range-dependent  sound  speed  profile  (SSP)  between  a 
fixed  source  and  an  array  of  receivers.  The  characteriza- 
tion of  the  SSP  from  acoustical  measurements  generally 
involves  inverse  techniques  in  order  to  infer  the  acoustic 
properties  of  the  medium  from  the  pressure  field  measured  at 
the  recei'v''ers. 

Due  to  the  complexity  of  the  ocean  the  inverse  problem 
is  most  often  non-linear  and  underdetermined.  Classicial 
acoustic  tomography  solves  this  problem  by  linearization. 
The  tomographic  method  is  able  to  estimate  the  perturbations 
of  the  sound  speed  field  by  comparing  the  measured  travel 
times  of  particular  rays  with  those  computed  numerically 
from  a  reference  sound  speed  field  and  a  raytrace  or  a 
normal  mode  algorithm  (Munk  and  Wunsch,  1978)  .  The 
procedure  provides  maps  of  the  perturbations  in  the  sound 
speed  field  and  indicates  how  different  the  actual  field  is 
from  the  one  used  as  a  reference.   If  a  large  enough  number 


of  ray  paths  are  used  in  the  computations,  the  spatial 
resolution  in  the  map  of  the  sound  speed  field  may  be  better 
than  the  one  obtained  from  discrete  CTD  measurements  (Howe, 
1986)  . 

Matched  field  processing  is  a  different  type  of  inverse 
method  which  was  first  proposed  as  a  method  to  locate  an 
acoustic  source  in  the  ocean.  The  principle  is  to  compare 
the  measured  complex  acoustic  pressures  at  a  vertical  array 
with  those  computed  from  an  acoustic  model  using  various 
positions  of  the  target  source  (Bucker,  1976) .  The 
procedure  generates  a  function  which  is  a  measure  of  the 
likelihood  between  the  actual  acoustic  pressure  field 
created  by  the  source  (unknown  position)  and  a  replica 
pressure  field  generated  from  an  estimate  of  the  source 
location. 

This  study  is  an  attempt  to  use  matched  field  processing 
as  an  alternate  tool  to  solve  the  inverse  problem  in 
acoustic  tomography.  Given  the  position  of  the  source  and 
the  receivers,  matched  field  processing  is  employed  to 
compare  the  true  complex  acoustic  pressures  at  the  receiving 
array  with  the  ones  computed  from  an  acoustic  propagation 
model  and  various  sound  speed  fields.  Computer  simulations 
are  used  to  demonstrate  the  performance  of  various 
estimators  under  different  signal-to-noise  ratios  and  noise 
correlation  matrix  structures. 


Chapter  II  provides  a  theoretical  presentation  of  the 
likelihood  estimators  which  are  used  in  this  study.  It  also 
illustrates  how  the  noise  is  modeled  and  added  to  the 
simulated  data.  The  simulations  are  shown  in  Chapter  III 
for  various  conditions  of  noise.  The  first  simulations  deal 
with  the  simple  case  where  only  one  parameter  of  the  medium 
is  unknown  (e.g.,  SOFAR  axis  depth,  surface  sound  speed, 
single  frontal  boundary)  and  where  the  noise  is  absent.  The 
next  cases  are  applied  to  the  localization  of  an  eddy  in  a 
noisy  medium.  Situations  of  both  spatially  uncorrelated 
noise  and  correlated  noise  were  examined.  Comparison  of  the 
estimators  is  provided  in  Chapter  IV;  the  spreading  of  each 
likelihood  function  about  the  true  value  is  examined  in  more 
detail  under  several  conditions  of  noise  power  and  noise 
correlation.  Also  discussed  is  the  problem  of  incomplete  or 
poor  knowledge  of  the  other  environmental  parameters,  e.g., 
incorrectly  specifying  the  bottom  absorption  property,  and 
how  this  may  introduce  inconsistency  in  the  procedure  and 
lead  to  a  decrease  in  the  efficiency  of  the  likelihood 
functions. 


II.   ACOUSTIC  TOMOGRAPHY  AND  MATCHED  FIELDS 

A.   PRINCIPLES  OF  MATCHED  FIELD  PROCESSING 

Classical  beamforming  for  plane  waves  is  obtained  by 
measuring  the  maximum  likelihood  between  the  actual  value  of 
the  complex  signal  at  each  hydrophone  and  the  values 
computed  from  an  expected  bearing  (Ziomek,  1985) . 

In  a  similar  fashion,  the  distance  to  a  target  in  the 
near  field  can  be  estimated  by  comparing  the  actual  values 
with  those  computed  for  different  distances.  The  range  is 
assumed  to  be  correct  when  both  sets  of  values  correspond  to 
the  same  wave  front  curvature  (Ziomek,  1985)  . 

Matched  field  processing  has  been  used  traditionally  to 
find  the  location  of  an  acoustic  source  in  a  shallow  water 
environment.  The  general  principle  is  to  store  the  values 
of  the  received  signal  (amplitudes  and  phase)  at  each 
element  of  the  array  and  then  compare  them  with  theoretical 
values  computed  for  different  possible  positions  of  the 
target  source.  The  true  location  of  the  emitter  is 
determined  when  both  fields  match  (Bucker,  1976;  Baggeroer, 
Kuperman  and  Schmidt,  1988). 

Different  criteria  may  be  used  to  measure  the  likelihood 
or  degree  of  matching.  Each  one  generates  a  different 
function  which  is  generally  well  adapted  for  a  particular 


type  of  noise.   Matched  field  detection  is  consistent  when 
the  medium  is  completely  determined. 

Matched  field  tomography  deals  with  the  inverse  problem. 
Given  that  the  source  location  is  known,  the  purpose  of  the 
procedure  is  to  estimate  the  medium  characteristics, 
particularly  the  range-dependent  sound  speed  profile.  In 
this  case,  the  replica  or  estimated  fields  are  built  from 
many  sound  speed  profiles  and  one  tries  to  match  them  with 
the  measured  one. 

B.   THEORY  IN  NOISE-FREE  CONDITIONS 
1 .   Bucker  Method 

According  to  Bucker  (1976),  the  following  "detection 
factor"  may  be  used  as  a  measure  of  the  difference  between 
the  exact  and  the  replica  fields: 

NR   NR 
y     "    <KS.,  •>  KR, 


_  j=l  k=I^j 


'jk   ^":k 


BUCK  =  -  "  ^  "^^ (2.1) 


where  terms  are  defined  as  follows: 


KR  =  spatial  autocorrelation  matrix  of  one  replica 
field, 

KS  =  spatial  autocorrelation  matrix  of  the  actual 
field, 

NR  =  number  of  hydrophones  in  the  array, 

F   =  scaling  factor  to  insure  a  result  between 
0  and  1, 


<>  =  time  average  used  to  remove  the  component  due 
to  the  noise, 

=  complex  conjugate 
The  matrices  KS  and  KR  are  defined  by 

KSjk  =  EjEk-  (2.2) 

KRjk  =  E'jE'k-  (2.3) 

where  p-;  and  p'j  denote,  respectively,  the  complex  envelope 
of  the  acoustic  pressure  and  the  replica  pressure  at 
hydrophone  j  ;  similarly  for  hydrophone  k.  As  demonstrated 
by  Bucker  (1976),  it  is  convenient  to  define  the  correlation 
matrices  from  the  complex  envelopes  of  the  signal  because 
the  rapidly-varying  time  component  is  removed  from  the 
computations.  These  complex  envelopes  are  easily  obtained 
by  processing  the  incoming  signal  through  a  classical 
quadrature  demodulator  (multiplication  by  a  sine  wave 
followed  by  a  low-pass  filter) . 

In  the  absence  of  noise,  the  time  average  is 
unnecessary  (KS  is  time  independent)  and  the  normalized 
detection  factor  becomes: 
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This  factor  is  similar  to  the  classical  correlation 
coefficient  of  two  random  variables.  The  expression  given 
above  does  not  use  the  diagonal  elements  of  the  matrices  and 
can  be  interpreted  as  the  output  of  a  regular  beamformer. 
Its  value  is  one  in  the  case  of  complete  equality  of  the 
fields  (KR  =  KS)  .  We  note  also  that  the  Bucker  detection 
factor  is  one  when  the  matrices  KR  and  KS  are  proportional, 
as  a  consequence  of  the  Schwarz  inequality. 
2 .   Heitmeyer  Method 

Heitmeyer  (1984)  defined  another  detection  factor 
which  he  called  the  "source  location  ambiguity  function." 
Its  value  is  given  by  the  following  expression: 
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and  may  be  rewritten  in  terms  of  phases  and  magnitudes: 
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It   can   easily   be   seen   that   the   function   is 
unnormalized.   From  the  inequality: 
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an  upper  limit  is  found: 
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The  ambiguity  function  is  always  bounded  by  a 
quantity  that  may  be  considered  as  the  average  power 
detected  at  each  hydrophone. 

When  the  replica  and  the  measured  fields  match 
exactly,  the  expression  reduces  to  the  following: 
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which  yields 
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The  inequality  demonstrated  previously  becomes  an  identity 
when  the  actual  and  replica  fields  are  identical. 


In  order  to  obtain  a  normalized  ambiguity  function, 
we  divide  Equation  (2.6)  by  this  upper  limit. 
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3 .   Relation  between  Bucker  and  Heitmeyer  Methods 

Using  the  definition  of  the  spatial  autocorrelation 
matrix,  Equation  (2.4)  may  be  written: 
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If  we  allow  the  subscripts  j  and  k  to  be  equal,  this 
expression  becomes: 
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hence. 
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or,  by  changing  the  subscripts, 
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which  is  exactly  the  normalized  ambiguity  function  defined 
by  Heitmeyer  in  Equation  (2.11). 

More  generally,  by  developing  the  expressions  of 
both  functions,  we  can  derive  the  following  relation  between 
the  Bucker  and  Heitmeyer  definitions: 
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This  relation  is  not  a  simple  proportionality  ratio 
because  it  changes  with  the  replica  fields  P'j- 
4 .   Center  of  Gravity  Method 

Other  likelihood  estimations  can  be  developed  using 
any   function  which  has   a  maximum   in  case   of  perfect 
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matching.   The  center  of  gravity  method  is  a  procedure  which 

solves  the  problem  from  a  more  mathematical  viewpoint. 

In  the  complex  plane,  vectors  P  and  P'  define  two 

sets   of  NR  points,   where  NR  represents  the  number  of 

receivers  of  the  array.   The  components  of  the  points  are 

the   real   and   imaginary   parts   of   the   complex   acoustic 

pressures: 

-v  , 

set  P  is  composed  of  points  (p-j_  cos  ii,Pi  sin  (fj_) 

->- 
set  P'  is  composed  of  points  (p'j_  cos  -*i/P'i  sin  -  '■[) 

In  order  to  compute  a  detection  factor,  we  calculate 
the  euclidean  distance  between  the  centers  of  gravity  G  and 
G'  of  both  sets  (Figure  2.1).  This  distance  is  inversely 
related  to  the  likelihood  of  the  fields. 

With  a  sum  of  weights  of  1,  the  points  G  and  G'  are 
given  by  their  coordinates: 
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(2.20) 


The  distance  becomes 


D=  ((Xa-Xc')2  +  (Y^-Ycj')2)  V2  (2.21) 
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We  then  normalize  the  quantity  in  order  to  obtain 
unity  for  complete  matching. 

DN  =  1  -  D/Dj^ax  (2.23) 

where  Dj^^^  i^  the  largest  unnormalized  distance  among  all 
the  replicas. 

DN  is  only  an  estimation  of  likelihood.  Although  it 
is  possible  for  two  different  sets  of  points  to  have  the 
same  center  of  gravity,  if  they  are  concentric,  this  does 
not  occur  in  the  simulations  and  the  method  keeps  its 
consistency.  We  will  see  later  that  this  distance  function 
may  lead  to  high  secondary  lobes  and  thus  is  not  always 
reliable. 

C.   THEORY  IN  PRESENCE  OF  NOISE 

Although  it  is  possible  to  use  the  former  expressions 
when  noise  is  present  which  contaminates  the  signal,  the 
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following  two  functions  are  more  specifically  suited  for  use 
in  the  presence  of  noise. 

Johnson  (1982)  previously  demonstrated  the  equivalence 
between  the  problem  of  bearing  determination  and  the 
estimation  of  the  spectrum  of  a  signal.  Due  to  this 
similarity,  all  modern  spectral  estimation  algorithms  apply 
equally  to  beamforming  and  matched  field  processing. 

Although  many  functions  may  be  used,  as  for  example, 
MUSIC  (Schmidt,  1981)  or  linear  predictor  (Johnson,  1982), 
we  will  particularly  emphasize  the  Bartlett  and  Maximum 
Likelihood  parameters  which  are  two  powerful  estimators  in 
target  location  problems.  These  estimators  are  especially 
useful  in  noisy  conditions,  because  Baggeroer  and  his 
colleagues  (1988)  showed  that  they  reduce  to  the  Bucker  or 
Heitmeyer  structures  when  the  noise  is  absent. 

1.   Bartlett  Method 

This  method  comes  directly  from  spectral  estimation 
theory.  Baggeroer,  Kuperman  and  Schmidt  (1988)  demonstrated 
that  the  power  output  of  a  Bartlett  beamformer  could  be 
written  in  the  following  quadratic  form: 

BART  =  W-  KT  W  (2.24) 

where  W  represents  the  normalized  velocity  potential  vector 
of  the  replica  field  and  KT  is  the  total  spatial  correlation 
matrix  of  the  signal  embedded  in  noise. 
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Due  to  the  proportionality  between  the  velocity 
potential  and  the  pressure  field,  the  equivalent  expression 
will  be  used: 

->■        ^  -V        _v 

BART  =(P'/|P'I)*KT(PV|P'I)  (2.25) 

where  P'  is  the  complex  acoustic  pressure  vector  of  the 
replica  at  the  array.  Under  the  condition  of  perfect 
matching  and  no  noise, 


NR    p    -  2 
BART  =yiD.  |=|P:  (2.26) 

1=1-^ 


which  is  the  summation  of  all  signal  powers  among  the 
hydrophones. 

In  order  to  normalize  the  function,  we  will  divide 
the  estimator  by  its  largest  value.  For  comparison  between 
the  different  methods,  we  will  focus  on  the  width  of  the 
main  lobe  rather  than  its  absolute  value. 

If  pj^  and  n-L  denote,  respectively,  the  signal  and 
the  noise  pressure  at  hydrophone  i,  the  spatial  correlation 
matrix  has  the  value: 


KT 


ij  =  E((pi+ni) (Ej+nj) •)  (2.27) 


where  E(   )  denotes  the  expectation  operation. 
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When  the  signal  and  the  noise  are  uncorrelated,  the 
matrix  reduces  to  the  simple  sum  of  signal  and  noise 
matrices: 

KTij  =  E(EiEj-)  +  E(ninj-)         .  (2.28) 

KTij  =  KSij  +  KNij  (2.29) 

These  matrices  are  hermitian  and  at  least 
semidef inite . 

2 .   Maximum  Likelihood  Method 

In  spectral  estimation,  the  Maximum  Likelihood 
method,  also  called  Capon's  method  or  the  minimum  variance 
method,  is  used  to  compte  the  power  spectral  density  of  a 
random  process  (Kay,  1988).   Its  expression  is  given  by: 

PMv(f)  =  (e"  R^x"^  e)-l  (2.30) 

where: 

f  =  frequency, 

^xx  ~  time  correlation  matrix  of  the  process, 

e  =  vector  whose  i^^  component  is  e^^^, 

H  =  transposition  of  the  conjugate  matrix. 

In  a  similar  way  the  output  of  a  Maximum  Likelihood 
beamformer  is  defined: 
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MLM  =  (W-  KT"1  W)"l  (2.31) 

where  W  and  KT  have  previously  been  defined. 

Following  the  procedure  of  Equation  (2.25),  Equation 
(2.31)  is  modified: 

MLM  =  ((P'VIP'I)  KT"1  (P'/|P'  I  )  )"■'■  (2.32) 

In  the  absence  of  noise,  KT  reduces  to  the  spatial 
autocorrelation  of  signal  only,  KS  (Equation  (2.29)).  As 
can  be  seen  when  the  array  is  composed  of  two  hydrophones, 
the  matrix  is  generally  singular  and  has  no  inverse.  The 
calculation  is  made  possible  by  adding  a  small  amount  of 
noise  to  the  diagonal.  As  with  the  Bartlett  estimator,  the 
Maximum  Likelihood  factor  will  be  divided  by  its  largest 
value  for  normalization. 
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III.   SIMULATION  OF  ACOUSTIC  TOMOGRAPHY 

Matched  field  processing  has  the  same  form  as  classical 
beamforming.  However,  instead  of  comparing  the  actual  field 
vector  with  a  plane  wave  replica  vector,  we  will  try  to 
match  the  actual  vector  with  a  vector  computed  from  an 
acoustic  propagation  model.  Measured  data  will  also  be 
simulated  with  the  same  code,  then  embedded  or  not  in  noise 
depending  on  the  scenario  under  investigation. 

A.   PROCEDURE 

1 .   Description  of  the  Simulation 

The  receiver  is  modeled  as  a  vertical  array  and  is 
assumed  to  be  composed  of  2  0  hydrophones  evenly  distributed 
between  550  m  and  1500  m.  The  source  is  located  100  km  from 
the  receiver  at  a  depth  of  1000  m.  It  emits  a  pure  sine 
wave  (tonal)  centered  at  100  Hz.  Due  to  the  inherent 
limitation  of  vertical  angles  in  the  parabolic  approximation 
(Ziomek,  1985) ,  the  transmitter  was  selected  to  have  a 
beamwidth  of  40°.  The  bottom  is  5000  m  deep  and  is  assumed 
to  be  flat  and  fully  absorbing.  This  choice  was  made  to 
speed  up  the  calculations  and  is  not  a  restrictive 
assumption.  It  assumes  all  the  energy  propagates  by  water- 
borne  paths. 
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2 .  Acoustic  Propagation  Model 

Several  models  could  have  been  used  to  simulate  the 
acoustic  fields,  for  example,  normal  mode  theory  or  the 
parabolic  equation  (PE) .  The  PE  model  was  used  because  it 
is  more  suitable  for  a  deep  water  simulation;  for  example, 
an  ocean  bottom  of  5000  m  allows  almost  670  propagating 
modes  at  100  Hz  and  would  have  been  computationally 
intensive  using  normal  mode  theory. 

The  classical  PE  approximation  with  split-step 
Fourier  transform  (Coppens,  1982)  was  available  in  the 
Environmental  Acoustic  Research  Group  package  of  models 
resident  at  NPS .  The  source  code  was  slightly  modified  to 
save  the  complex  acoustic  pressures  at  the  hydrophones  in  a 
file.  The  measured  and  the  replica  complex  pressure  fields 
were  then  stored  in  order  to  run  the  simulation  programs. 

3 .  Simulation  of  Noise 

Noise  was  added  to  the  measured  data  in  order  to 
produce  a  realistic  problem  and  to  study  the  behavior  of  the 
estimators  in  different  environments.  Following  the 
procedure  described  by  Porter,  Dicus  and  Fizell  (1987), 
noise  was  introduced  by  the  mean  of  its  spatial  correlation 
matrix.  This  procedure  is  better  than  just  altering  the 
data  with  random  noise  because  the  probability  density 
function  of  the  noise  is  difficult  to  estimate.  Moreover, 
all  the  estimators  considered  were  written  in  terms  of 
correlation  matrices. 
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Ambient  noise  falls  in  two  categories: 

-  uncorrelated  noise.  Its  correlation  matrix  is 
proportional  to  the  identity  matrix  and  the 
proportionality  factor  is  an  indicator  of  the  noise 
power. 

-  correlated  noise.  In  this  case  the  matrix  has  non  zero 
terms  outside  the  diagonal,  but  is  nevertheless  an 
hermitian  matrix. 

Several  attempts  to  measure  the  coherence  of  ambient 
noise  in  the  ocean  have  been  conducted  during  the  past 
years.  One  of  them  was  made  from  the  Trident  Vertical  Array 
and  is  described  by  Urick  (1984).  Figure  3.1  depicts  the 
results  of  this  study  and  has  been  used  to  generate  a  model 
of  the  noise  correlation  matrix. 

The  matrix  was  modeled  in  the  following  way: 

KNij  =  r2  g-  i  i-j  I  (3.1) 

where  z^  depends  on  the  noise  power  and  ■  is  a  factor  which 
indicates  how  fast  the  coherence  falls  off  outside  the 
diagonal.  The  larger  ;  is,  the  more  uncorrelated  the  noise 
is,  i.e.,  its  spatial  correlation  scale  becomes  shorter. 

Due  to  the  spacing  between  the  receivers  in  the 
array  and  the  frequency  (100  Hz)  used  in  the  simulation, 
Figure  3.1  shows  that  =  =  1.7  is  consistent  with  the 
observed  ambient  noise  correlation. 

Although  the  spatial  correlation  matrix  of  the 
noise,  KN,  is  generally  a  complex  hermitian  matrix,  the 
analytic  modeling  shown  in  Equation  (3.1)  describes  a  real 


20 


symmetric  matrix.  As  explained  by  Cox  (1973),  this 
approximation  is  valid  in  the  special  case  of  zero  time 
delay,  i.e.,  when  it  is  assumed  that  the  noise  is  in  phase 
among  all  the  hydrophones  of  the  vertical  array.  This 
assumption  is  relatively  consistent  for  low  frequency  noise. 
Below  150  Hz  the  noise  is  principally  due  to  distant 
shipping  and  arrives  mainly  from  the  horizontal. 

B.   ONE  DIMENSIONAL  PROBLEM  (NOISE-FREE  CONDITIONS) 

In  the  following  simulation,  the  shape  of  the  sound 
speed  profile  is  the  only  unknown.  If  the  profile  is 
digitized  in  1  m  intervals,  then  for  a  water  depth  of  5000 
m,  one  would  have  to  determine  5000  values  to  match  the 
complete  sound  speed  profile.  Such  a  procedure  would  lead 
to  an  unmanageable  number  of  computations,  especially  if  one 
tries  to  match  a  large  number  of  replica  fields.  However, 
as  we  are  only  interested  in  demonstrating  the  feasibility 
of  the  procedure,  we  will  begin  by  studying  the  simple  cases 
where  only  one  or  two  points  of  the  sound  speed  profile  are 
unknown . 

1 .   Determination  of  a  SOFAR  Axis  Depth 

We  initially  start  with  a  bi-gradient  sound  speed 
profile  having  a  sound  speed  minimum  at  1000  m  (Figure  3.2). 
Two  replica  profiles  are  considered  wherein  the  SOFAR  axis 
depth  is  altered  by  +/-  200  m  (step  10  m)  .  Note  that 
because  the  surface  and  bottom  sound  speeds  are  unchanged, 
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the  change   in  axial  depth  results   in  a  change  in  the 
gradients  of  both  the  upper  and  lower  segments  of  the  SSP. 

Three  estimation  techniques,  the  Bucker,  Heitmeyer, 
and  center  of  gravity  methods  are  utilized  to  determine  the 
true  depth  of  the  sound  speed  minimum.  These  estimators  are 
computed  from  Equations  (2.4),  (2.11)  and  (2.23).  The 
estimated  depth  of  the  SOFAR  axis  is  found  when  an  estimator 
shows  a  peak  with  a  detection  factor  of  1.  The  results  are 
shown  in  Figures  3.3,  3.4  and  3.5;  each  estimator 
demonstrates  a  different  behavior. 

The  Bucker  detection  factor  and  the  Heitmeyer 
ambiguity  function  both  indicate  a  maximum  at  the  true 
location  of  1000  m.  However  a  strong  side  lobe,  centered  at 
1040  m,  indicates  these  two  detection  factors  are  not  robust 
enough  to  provide  an  unambiguous  selection  of  the  SOFAR  axis 
depth.  In  addition,  strong  secondary  side  lobes  are  also 
present. 

The  center  of  gravity  method  (Figure  3.5)  proves  to 
be  a  better  estimator  for  this  situation.  The  main  lobe  is 
much  narrower  and  no  other  lobes  exist.  For  a  noise-free 
ocean,  this  is  the  best  estimator  among  the  three  to  solve 
this  particular  problem. 

2 .   Determination  of  Surface  Sound  Speed 

In  order  to  mimic  typical  seasonal  or  spatial 
changes  in  the  SSP,  the  surface  sound  speed  was  permitted  to 
fluctuate  by  +/-  5  m/s  about  a  mean  value  (step  1  m/s) .   For 
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this  situation  the  SOFAR  axial  depth  was  fixed  at  1000  m. 
Hence,  only  the  upper  gradient  changes  as  seen  in  Figure 
3.6.  We  seek  an  estimator  that  will  match  the  true  surface 
sound  speed  (solid  line  in  Figure  3.6).  Using  the  same 
three  estimation  techniques  as  above,  the  results  are  shown 
in  Figures  3.7,  3.8  and  3.9.  For  this  situation,  we  observe 
a  nearly  identical  behavior  of  the  Bucker  and  the  Heitmeyer 
functions,  both  of  which  have  moderate  side  lobes  at  about 
0.7.  As  before,  the  center  of  gravity  estimator  remains  the 
best  without  any  ambiguity  due  to  the  presence  of  secondary 
lobes. 

3 .   Determination  of  an  Acoustic  Frontal  Boundary 

To  model  the  presence  of  an  acoustic  front  two 
different  sound  speed  profiles  are  introduced,  one  50  km 
from  the  other.  The  parabolic  equation  model  was  utilized 
in  this  range-dependent  problem  with  the  position  of  the 
front  (i.e.,  the  range  at  which  the  second  SSP  is 
encountered)  allowed  to  vary  by  +/-  2  0  km  about  the  true 
position  (step  1  km).  Figure  3.10  provides  an  illustration 
of  the  SSP  setup.  The  Bucker  and  the  Heitmeyer  methods 
correctly  solve  this  problem  with  relatively  narrow  main 
lobes,  as  shown  in  Figures  3.11  and  3.12.  However,  the 
center  of  gravity  method  does  not  perform  as  well  (Figure 
3.13).  Although  it  is  able  to  locate  the  correct  value, 
significant  sidelobes  are  present  which  could  lead  to  an 
ambiguity  if  too  small  a  detection  threshold  were  chosen. 
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This  last  procedure  is  obviously  not  suitable  for  these 
conditions. 

4 .   Comments  on  the  Estimators 

In  order  to  completely  appreciate  the  consistency  of 
each  of  the  previous  estimators  for  the  case  where  only  one 
factor  is  unknown,  it  is  important  to  study  their  response 
to  a  variety  of  configurations. 

a.   Influence  of  the  Number  of  Hydrophones 

In  the  problem  of  determining  the  depth  of  the 
SOFAR  axis,  the  array  was  composed  of  2  0  hydrophones.  It  is 
possible  to  run  the  same  simulation  by  using  only  a  fraction 
of  the  receivers.  Plots  of  the  Bucker  estimator  for  four 
different  numbers  of  hydrophones  are  shown  in  Figure  3.14. 
Although  the  difference  is  small  when  the  number  of 
hydrophones  is  reduced  from  20  to  five,  the  output  of  an 
array  composed  only  of  two  receivers  changes  drastically. 
When  the  number  of  hydrophones  is  this  small,  the  side  lobes 
may  have  an  amplitude  of  the  same  order  as  the  main  lobe,  a 
situation  which  leads  to  full  ambiguity.  The  number  of 
hydrophones  is  thus  an  important  parameter  which  must  always 
be  more  than  some  minimum  value.  This  value  is 
unfortunately  dependent  upon  the  actual  problem  and  the 
depth  of  the  array  relative  to  the  axial  depth.  Receivers 
which  do  not  intercept  much  of  the  acoustic  energy  can  be 
easily  omitted  but  those  which  contain  significant  amplitude 
and  phase  information  should  be  retained. 
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b.  Influence  of  the  Frequency 

Typically  as  a  result  of  beamforming,  the  main 
lobe  becomes  narrower  as  the  frequency  of  the  array 
increases.  The  same  phenomenon  can  be  observed  in  Figure 
3.15,  wherein  the  width  of  the  estimator  peak  is  also  a 
function  of  the  frequency.  However,  a  trade-off  exists 
between  the  desired  resolution  (width  of  the  detection  peak) 
and  the  computation  time  of  the  PE  model  which  increases 
rapidly  with  frequency.  Also  if  higher  frequencies  (kilo- 
Hertz  range)  were  used,  the  signal  would  be  limited  by 
absorption  which  would  result  in  lower  signal-to-noise 
ratios. 

c.  Importance  of  Array  Position 

The  depth  of  the  array  is  of  minimal  importance 
when  working  in  shallow  water  because  the  entire  water 
column  is  nearly  insonified.  Such  is  not  the  case  in  deep 
water  where  shadow  zones  exist  with  relatively  low  signal 
levels.  Based  on  the  depth  of  the  source,  a  first  guess  of 
the  depth  to  position  the  array  would  be  to  place  it  where 
the  signal  may  be  expected  to  occur  with  a  high  level,  for 
example,  in  the  vicinity  of  the  SOFAR  axis  or  near  a 
convergence  zone.  The  choice  will  obviously  be  dependent  on 
the  profile  shape.  " 
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C.   TWO  PARAMETER  PROBLEM  (NOISY  CONDITIONS) 
1.   Localization  of  an  Eddy 

The  previous  section  has  shown  that  the  estimators 
are  generally  able  to  find  the  correct  value  in  the  case  of 
a  simple  unknown  and  a  noise-free  medium.  The  same  kind  of 
simulation  may  be  run  when  two  parameters  are  to  be 
determined.  To  test  the  ability  of  the  various  estimators 
to  deal  with  a  two  dimensional  problem,  we  will  examine 
their  ability  to  locate  an  eddy  assumed  to  be  present 
between  a  source  and  a  receiving  array.  The  sound  speed 
profiles  inside  and  outside  the  eddy  are  known.  Thus  the 
only  unknowns  are  the  borders  of  this  perturbation  of  the 
sound  field.  An  eddy  20  km  in  diameter  is  positioned  40  km 
to  60  km  from  the  source.  The  replica  fields  are  computed 
by  scanning  the  limits  from  35  km  and  45  km  for  the  border 
closest  to  the  source,  and  from  55  km  to  65  km  for  the 
farther  boundary.  Replications  at  1  km  interval  were  made. 
We  will  thus  try  to  match  the  simulated  measured  data  with 
121  replica  fields.  Figure  3.16  presents  the  true  location 
of  the  eddy  in  this  simulation. 

Figure  3.17  shows  the  disposition  of  source  and 
array  with  regard  to  the  energy  field  for  the  true  location 
of  the  eddy.  The  array  lies  on  the  SOFAR  axis,  almost  3  0  km 
beyond  a  convergence  zone.  From  this  plot,  we  can  expect  a 
high  level  of  signal  from  the  channel  propagation  and  large 
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differences   in   phase   due   to  multipath   propagation   (RR 
acoustic  rays) . 

Before  processing,  the  measured  data  are  imbedded  in 
noise.  This  noise  is  introduced  by  the  spatial  autocorrela- 
tion matrix  described  earlier.  The  principal  assumption  of 
this  simulation  is  that  all  correlation  matrices  are 
completely  known.  For  an  actual  situation,  this  may  not  be 
true  but  it  is  still  possible  to  estimate  the  total  matrix 
of  noise  from  the  set  of  measured  data. 

Because  of  its  close  similarity  to  the  Bucker 
detection  factor,  the  Heitmeyer  function  will  be  omitted 
from  further  analysis.  The  Bartlett  and  the  Maximum 
Likelihood  functions  are  introduced  for  these  simulations, 
and  it  will  later  be  seen  that  these  two  estimators  are  well 
suited  for  conditions  where  the  signal-to-noise  ratio  is 
low. 

2 .   Siqnal-to-Noise  Ratio 

The  signal-to-noise  ratio  is  a  parameter  which 
depends  on  the  relative  powers  of  the  signal  and  the  noise 
at  the  array.  Following  the  procedure  of  Ziomek  (1985),  the 
signal-to-noise  ratio  can  be  written  in  terms  of  the  noise- 
free  signal  and  the  spatial  autocorrelation  matrix  of  the 
noise: 

NR     NR 

.1     i  HiPj- 
SNR   =    20    Logio  ^g-^^ (3.2) 

)         ^    KN. . 
i=l   j=l      "3 
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The  numerator  of  this  expression  may  be  interpreted 
as  the  summation  of  all  the  elements  of  the  signal  spatial 
correlation  matrix.  Similarly,  the  denominator  represents 
the  sum  of  the  entries  of  the  noise  correlation  matrix. 

The  SNR  was  computed  for  several  values  of  the  noise 

power,  o^ ,  and  the  correlation  parameter,  :<,     that  generate 

new  values  of  KNj^j   (Equation   (3.1)).    The  results  are 

presented  in  Table  3.1.   The  expression  above  shows  that  the 

SNR   is   reduced   when   the   denominator   of   the   argument 

NR  NR 
increases,   i.e.,  when  the  double  summation        ^^^ii  i^ 

i=l  j=l 
large  due  to  a  significant  increase  in  the  power  of  the 

noise,  -  2,  or  a  slow  decay  in  the  correlation  between  the 

hydrophones,  B.    By  using  the  results  presented  in  Figure 

3.1  and  the  analytic  expression  of  the  noise  matrix  given  in 

Equation   (3.1),   the   correlation   of   the   noise   will   be 

considered  high  when  the  factor  ,:  is  less  than  0.57  and  low 

when  it  exceeds  2.0. 

TABLE  3 . 1 
SIGNAL  TO  NOISE  RATIO  AS  A  FUNCTION  OF  ^  ^  p^^^^    3 


o2        ' 

0.! 

57 

1.0 

1.7 

2.0 

2.2 

10-19 

13 

dB 

16 

dB 

20  dB 

21 

dB 

21  dB 

10-19 

-1 

dB 

3 

dB 

6  dB 

7 

dB 

7  dB 

10-18 

-7 

dB 

-3 

dB 

0  dB 

1 

dB 

1  dB 

10-18 

-21 

dB 

-17 

dB 

-14  dB 

-13 

dB 

-12  dB 

10-1-7 

-27 

dB 

-23 

dB 

-20  dB 

-19 

dB 

-19  dB 

10-16 

-47 

dB 

-43 

dB 

-39  dB 

-39 

dB 

-38  dB 

10-18 

-67 

dB 

-63 

dB 

-60  dB 

-59 

dB 

-58  dB 
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In  the  following  sections,  the  performance  of 
various  estimators  will  be  examined  for  several  conditions, 
including  different  cases  of  noise  power  and  noise 
correlation. 

3 .   Bucker  Detection  Factor 

A  simulation  using  the  Bucker  detection  factor  was 
run  for  the  case  of  eddy  localization  under  both  noise-free 
and  noisy  conditions. 

a.   Noise-free  Conditions 

The  simulation  was  run  by  setting  the  power  of 
the  noise,  :^,  to  zero.  The  autocorrelation  matrix  of  the 
noise  is  then  just  the  null  matrix  and  the  total  matrix 
reduces  to  one  of  signal  only,  as  indicated  by  Equations 
(3.1)  and  (2.29).  A  three  dimensional  plot  and  a  contour 
plot  of  the  detection  factor  are  shown  in  Figure  3.18.  The 
estimator  is  represented  by  a  surface  which  has  only  a 
single  maximum  positioned  at  the  correct  location  of  the 
eddy.  This  suggests  that  the  Bucker  method  is  able  to 
determine  the  true  location  of  the  eddy  in  a  noise-free 
environment. 

An  interesting  feature  of  the  plot  is  the 
symmetry  that  exists  around  both  diagonals  of  the  contour. 
Moving  on  the  principal  diagonal,  along  the  line: 

Y  =  X  +  20  (3.3) 
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is  equivalent  to  displacing  an  eddy  with  a  constant  diameter 
of  20  km.  The  small  intervals  between  the  contour  lines 
along  this  path  indicate  that  the  location  of  the  eddy  can 
be  determined  with  reasonable  accuracy,  once  we  know  its 
diameter. 

b.   Case  of  Uncorrelated  Noise 

Uncorrelated  noise  is  generated  by  choosing  a 
large  value  of  :  .  In  this  case,  the  correlation  falls  off 
rapidly  on  either  side  of  the  noise  matrix  diagonal.  For  a 
large  enough  : ,  the  noise  field  at  one  hydrophone  is 
completely  dissimilar  to  that  at  another  hydrophone  and  the 
noise  matrix  becomes  diagonal.  Simulations  were  run  with 
iS=  10  and  different  values  of  r^.  All  yielded  the  same 
results  as  in  Figure  3.19,  which  is  seen  to  be  identical  to 
Figure  3.18.  This  similarity  may  be  explained  by  recalling 
the  definition  of  the  Bucker  detection  factor  when  noise  is 
present: 


NR   NR 
y    y   KT.,   KR., 

BUCK  =  3-1  K-ij^j (2.4) 

NR   NR  ,  ,^  NR   NR  -,  /^ 

(  y     y   KR    KR   )l/2(  I  I     KT   KT   O''^' 

j=l  k=l7^j   ^^     3^      j=l  k=l?^j   ^^    ^^ 


where  the  total  correlation  matrix  is  given  by 


KTj)^  =  KSj}^  +  KNjk  (2.29) 
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or, 


KTjk  =  KSjk  +  -2  exp(-,|j-k|)  (3.4) 

As  can  be  seen,  this  expression  only  uses  the 
cross  terms  of  the  matrix.  When  6  is  large  enough,  the 
second  term  in  the  right  hand  side  of  Equation  (2.29)  is 
almost  negligible  for  every  value  of  the  noise  power,  ^2; 
thus  the  cross  terms  of  the  total  matrix  KT  reduce  to  the 
cross  terms  of  the  correlation  matrix  of  the  signal  KS . 

KTjk   ;   KSjk,   j  ^  k  (3.5) 

By  changing  the  value  of  -^ ,  the  power  of  the 
noise  is  modified,  but  the  new  diagonal  terms  do  not  play  a 
role  in  the  calculations.  The  Bucker  detection  factor  is 
thus  insensitive  to  perfectly  uncorrelated  noise;  in  this 
case,  the  performance  is  exactly  identical  to  the  one  in 
noise-free  conditions. 

c.   Case  of  Uncorrelated  Noise 

Any  combination  of  noise  power,  -2^  ^nd  noise 
correlation,  -,  yields  a  different  pattern  of  the  detection 
factor.  In  cases  for  which  the  spatial  correlation  of  the 
noise  is  high,  the  Bucker  method  may  still  be  able  to 
maximize  the  detection  factor  at  the  correct  location,  but 
the   absolute   value   of   the   peak  will   decrease   as   the 
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ambiguity  surface  becomes  flatter.   Figure  3.20  illustrates 
this  type  of  behavior. 

Whenever  both  '  and  "^  generate  a  low  signal-to- 
noise  ratio  (large  power,  z^  ,  or  small  correlation 
parameter,  .:  )  ,  the  procedure  fails  and  the  localization  of 
the  eddy  becomes  impossible  (see  Figure  3.21).  We  will 
quantify  the  effects  of  "2  and  i.  on  localization  below. 
4 .   Bartlett  Estimator 

The  same  simulations  as  above  were  run  using  the 
Bartlett  estimator  under  noise-free  and  noisy  conditions. 

a.  Noise-free  Medium 

In  a  generic  noise-free  environment,  the 
performance  of  the  Bartlett  estimator  is  similar  to  the 
Bucker  detection  factor,  as  shown  in  Figure  3.22.  As 
expected,  the  same  symmetry  along  the  diagonals  is  still 
present. 

b.  Correlated  Noise 

The  performance  of  the  Bartlett  estimator 
changes  significantly  as  -^  ^nd  ■  vary.  Figure  3.23 
provides  an  example  of  the  plot  for  r^  =  lo"^^  and  ^  =  1.7, 
where  the  true  location  is  found.  Figures  3.24  and  3.25 
illustrate  failures  of  the  method  due  to  a  weak  signal-to- 
noise  ratio  brought  about  by  strong  noise  and  highly 
correlated  noise,  respectively.  It  will  later  be 
established  that  the  usual  characteristics  of  the  noise  in  a 
deep  ocean  do  not  generally  lead  to  this  kind  of  ambiguity. 
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5 .   Maximum  Likelihood  Method 

The  MLM  estimator  was  calculated  using  the  same 
noise  conditions  as  for  the  Bartlett  function.  An 
examination  of  noise-free  conditions  is  not  possible  because 
of  the  singularity  of  the  correlation  matrix.  Equation 
(2.31)  shows  that  the  expression  of  the  Maximum  Likelihood 
estimator  requires  the   calculation   of   the   inverse  matrix 

MLM  =  (W-  KT"1  W)"l  (2.31) 

When  noise  is  absent,  the  matrix  KT  reduces  to  the 
correlation  matrix  of  the  signal  KS  which  is  generally 
singular. 

a.  Slightly  Correlated  Noise 

When  -2  =  10"-^^  and  ■  =  1.7,  the  method  gives 
better  results  than  the  Bartlett  estimator  with  relatively 
low  side  lobes  (compare  Figure  3.26  with  Figure  3.23). 

b.  Strongly  Correlated  Noise 

With  a  more  highly  correlated  ambient  noise 
(3  =  0.17),  as  depicted  in  Figure  3.27,  it  is  still  possible 
to  obtain  a  correct  location  of  the  eddy.  By  comparing  this 
plot  with  Figure  3.25,  we  note  that  the  Bartlett  estimator 
was  unsuccessful  in  this  case.  Nevertheless,  even  for  the 
MLM  technique,  a  very  low  signal-to-noise  ratio  will  result 
in  a  failure  as  shown  in  Figure  3.28. 
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Figure  3.2   Actual  and  Extreme  Replicas  of  the  SSP  used 
in  the  Simulation  to  Determine  SOFAR  Axis 
Depth 
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Figure    3 . 4 


SOFAR  Axis  Depth  Determination.   Heitmeyer 
Ambiguity  Function 
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Figure    3 . 5 
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Figure  3.6   Actual  and  Extreme  Replicas  of  the  SSP  used 
in  the  Simulation  to  Determine  the  Surface 
Sound  Speed 
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Figure  3.7   Surface  Sound  Speed  Determination.   Bucker 
Detection  Factor 
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Figure  3.8   Surface  Sound  Speed  Determination 
Heitmeyer  Ambiguity  Function 
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Figure  3.9   Surface  Sound  Speed  Determination.   Center 
of  Gravity  Method 
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Figure  3.10   Simulation  to  Determine  the  True  Location 
of  an  Acoustic  Front 
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Figure  3.11   Acoustic  Front  Determination 
Detection  Factor 
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Figure  3.12   Acoustic  Front  Determination.   Heitmeyer 
Ambiguity  Function 
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Figure  3.13   Acoustic  Front  Determination, 
of  Gravity  Method 


Center 


46 


ft: 
o 

rn    • 


UJ 

o 
m 


70 
10 

3 

2 


liydropltoting  !  ■ 
liyclropliOMR^  I  ■ 
hydropl  innf?  s  t 
liydr  oplionc5  ! 


800.0 


900.0  1000.0 

sornR  nxis  ucnii 


— t  •  ■ 

1100. 0 


1200.0 


Figure  3.14   SOFAR  Axis  Depth  Determination  for  Various 
Numbers  of  Hydrophones  in  the  Array 
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Figure  :.L7   Surface  Sound  Speed  Determination. 
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Figure  3.18   Eddy  Localization.   Bucker  Detection 
Factor.   Noise-Free  Conditions 
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Figure  3.20 


Eddy  Localization.   Bucker  Detection  Factor 
Illustrating  Condition  of  Highly  Correlated 
Noise  and  Moderate  Noise  Power,  '^  =  10"-^^, 
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Figure  3.22   Eddy  Localization.   Bartlett  Estimator, 
Noise  Free  Conditions 


55 


:5:"i.o  3fi.n 


37.0     30.0     n.n     in.u     -ii.o     42.0     13. 
CLU5E  LI  11  IT   []r  EDUY    (Kri) 


1!.n   15.0 


Figure  3.23 


Eddy  Localization.   Bartlett  Estimator 
Illustrating  Condition  of  Moderate  Noise 
Power  and  Moderate  Noise  Correlation, 
o2  =  10"16,  B  =  1.7 
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Figure  3.24  Eddy  Localization.  Bartlett  Estimator 
Illustrating  Condition  of  Strong  Noise 
Power  and  Moderate  Noise  Correlation, 
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Figure  3.25   Eddy  Localization.   Bartlett  Estimator 

Illustrating  Condition  of  Moderate  Noise 
Power  and  Strong  Noise  Correlation, 
o2  =  10~1^,  B  =  0.17 
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Figure  3.26   Eddy  Localization.   MLM  Estimator  Illus- 
trating Condition  of  Moderate  Noise  Power 


and  Moderate  Noise  Correlation, 

e=  1.7 


.2  _ 


10 


-16 


59 


■J5.0  36.0      37. n      3n.O      3D.0      10.0      41.0      12.0      13.0 

CLOSE   LIMIT  or   LDUY    (KM) 


11.0   15.0 


Figure  3.27   Eddy  Localization.   MLM  Estimator  Illus- 
trating Condition  of  Moderate  Noise 
Power  and  Strong  Noise  Correlation, 
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Figure  3.28 


Eddy  Localization.   MLH  Estimator  Illus- 
trating Condition  of  Strong  Noise  Power 
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IV.   ANALYSIS  OF  THE  SIMULATION 

From  the  previous  simulations,  one  sees  that  the 
performance  of  each  of  the  estimators  varies  significantly 
under  varying  signal-to-noise  ratios  or  source/receiver 
geometries  or  sound  speed  variations.  For  example,  the 
center  of  gravity  method  was  shown  to  be  the  best  in 
locating  the  SOFAR  axis  and  determining  the  surface  sound 
speed.  In  contrast,  this  procedure  was  the  least  successful 
in  the  acoustic  front  localization  problem.  Therefore  the 
efficiency  of  an  estimator  does  not  depend  only  on  the  type 
of  ambient  noise  but  also  on  the  particular  problem  being 
solved.  In  order  to  continue  focusing  on  a  realistic 
problem,  the  eddy  localization  problem  will  be  studied  in 
greater  detail.  The  condition  of  mismatching  will  be 
treated  separately. 

A.   COMPARATIVE  STUDY  OF  THE  ESTIMATORS 
1 .   Criterion 

To  facilitate  comparison  of  the  performance  of  the 
different  methods,  we  will  calculate  the  joint  central 
moments  of  the  different  estimators.  The  joint  central 
moment  provides  information  on  the  spread  of  the  detection 
factor  about  the  mean.  The  smaller  the  moment,  the  better 
the  estimate  of  the  parameter  of  the  ocean. 
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Using  the  definition  of  the  central  moment  of  a 
multiple  random  variable  as  defined  by  Peebles  (1987) : 


^nk  =  ^^  I  (x-X)^    (y-Y)^  fxyi^'Y)    dx  dy  (4.1) 


where: 

X,  Y  =  means  of  the  random  variables  X  and  Y, 

f^y   =  joint  probability  density  function  of  X  and  Y, 

n,  k  =  orders  of  the  central  moment. 

The  spreading  factor  for  the  MLM  method  is  defined 
as  the  second  order  central  moment  of  the  function  MLM(x,y): 


S  =    '  '  ^  MLM(x,y)   (x-40)2  (y-60)2  dx  dy  (4.2) 

'  X  V 


where  x  and  y  are  the  boundaries  of  the  eddy  we  are  looking 
for. 

Since  we  are  dealing  with  a  discrete  search  among 
parameter  values,  this  pseudo  variance  has  the  form: 


NR  NR 
S  =   y    y  MLMi-j  (44  +  1-40)2  (54  +  j-60)2         (4.3) 

i=l  j=l 


where  NF  denotes  the  number  of  possible  values  of  i  and  j . 

Equation  (4.3)  indicates  how  the  estimator  spreads 
around  the  true  frontal  boundary  values  x  =  40  and  y  =  60. 
It  is  a  global  measure  of  the  estimator  performance,  not 
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just  a  measure  of  the  main  lobe  width.  A  large  value  is 
associated  with  significant  spreading  which  indicates  a  poor 
performance,  even  if  the  true  position  is  actually  found. 
The  value  of  the  spreading  factor  will  be  also  large  when 
the  ambiguity  surface  has  a  large  mean  value  and  a  small 
amplitude  (case  of  significant  noise) .  Analogous  spreading 
factors  may  be  defined  for  the  Bucker  and  the  Bartlett 
functions. 

2 .   Efficiency  of  the  Estimators 

The  spreading  factor  defined  above  has  been  computed 
for  the  different  combinations  of  -^  ^nd  shown  in  Table 
3.1.  The  results  are  presented  in  Table  4.1  for  the  Bucker, 
Bartlett  and  MLM  methods.  From  this  table  it  is  possible  to 
compare  the  efficiency  of  each  technique  in  a  variety  of 
environments.  The  values  of  the  correlation  factor  .:  have 
been  chosen  to  stay  consistent  with  deep  water  measures. 
The  range  of  noise  powers,  -^  produced  signal-to-noise 
ratios  between  -67  dB  and  +21  dB. 

In  order  to  be  consistent  in  comparing  the  spreading 
factor  S  of  the  three  schemes,  it  is  convenient  to  modify 
the  expression  for  the  Bucker  detection  factor  as  defined  by 
Equation  (2.4)  by  dividing  it  by  its  largest  value  among  all 
the  replicas.  The  maximum  of  this  new  function  will  always 
be  one  in  all  situations  of  noise,  as  the  Bartlett  and  the 
MLM  estimators. 
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TABLE  4 . 1 

SPREADING  FACTOR  S  OF  THE  BUCKER  (NORMALIZED) ,  BARTLETT 
AND  MLM  ESTIMATORS  AS  A  FUNCTION  OF  :^    AND  - 

o    2      0.57         1.00         1.70         2.00         2.20 

SNR=+13dB  SNR=+16dB  SNR=+20dB  SNR=+21dB  SNR=+21dB 

10"!^   19,000  18,973  18,950      18,944  18,941 

25,924  25,902  25,884       25,880  25,877 

114  139  142          140  139 

SNR=-ldB  SNR=+3dB  SNR=+6dB  SNR=+7dB  SNR=+7dB 

5xl0"19   19,280  19,147  19,031  19,004  19,090 

26,471  26,364  26,272  26,250  26,239 

567  654  708  698  691 

SNR=-7dB  SNR=-3dB  SNR=OdB  SNR=+ldB  SNR=+ldB 

10"^^   19,632  19,364  19,133  19,078  19,050 

27,149  26,935  26,754  26,710  26,688 

1,131  1,382  1,411  1,391  1,377 

SNR=-21dB  SNR=-17dB  SNR=-14dB  SNR=-13DB  SNR=-12dB 

5>10"1S   22,400      21,081  19,946  19,671  19,532 

32,328      31,320  30,463  30,258  30,153 

5,499        6,657  6,766  6,670  6,603 

SNR=-27dB  SNR=-23dB  SNR=-20dB  SNR=-19dB  SNR=-19dB 

10"!^   25,773  23,188      20,950      20,407  20,131 

38,256  36,371      34,771       34,386  34,190 

10,630  12,732       12,883       12,687  12,570 

SNR=-47dB   SNR=-43dB  SNR=-39dB  SNR=-39dB  SNR=-38dB 

10"^^     "          54,747       37,381  32,723  30,280 

85,126      78,315  76,655  75,805 

66,998      72,147      70,060  69,050  68,404 

SNR=-67dB   SNR=-63dB   SNR=-60dB   SNR=-59dB   SNR=-58dB 

10~15nc  CO  00  cc 

00  «>  °^  °°  124,781 

oo  00  00  oo  124,584 

Note  1:   The  signal-to-noise  ratio  is  indicated  for  each 
entry  of  the  table,  followed  by  a  triplet  of  numbers  which 
represent,  respectively,  the  Bucker  (normalized) ,  Bartlett 
and  MLM  spreading  factors. 

Note  2:   The  value  =^  means  that  the  estimator  is  unable  to 
detect  the  true  location  of  the  eddy. 


65 


Table  4.1  provides  a  good  illustration  of  the 
performance  of  the  estimators  under  several  conditions  of 
noise  power  and  correlation.  Figures  4.1,  4.2  and  4.3  show 
logarithmic  plots  of  the  spreading  factor  S  versus  the  noise 
power  '^  ^ ;  each  curve  represents  a  value  of  the  correlation 
parameter  B .  It  is  thus  possible  to  determine  the  value  of 
S  for  any  pair  of  o^  and  z  .  As  is  obvious  from  Table  4.1, 
the  comparative  performance  of  each  method  is  mostly  a 
function  of  the  signal-to-noise  ratio  of  the  measure. 

For  SNR  less  than  -50  dB  all  methods  fail  to  detect 
the  true  location  of  the  eddy.  Nevertheless,  we  observe  the 
case  j^  =  10"^^  and  ■=  2.2,  where  the  Bartlett  and  MLM 
estimators  indicate  two  different  maxima  at  one.  Even  in 
this  case,  the  amplitude  of  the  functions  is  so  small  that 
the  spreading  factor  S  is  very  large. 

When  the  SNR  is  about  -40  dB,  the  Bucker  method  is 
the  most  efficient  estimator,  albeit  a  weak  one,  as  the 
spreading  remains  significant.  The  superiority  of  the 
scheme  increases  moreover  when  the  correlation  of  the  noise 
decreases. 

For  other  SNR  and  correlation  values,  the  MLM  method 
is  generally  the  most  efficient.  When  the  SNR  =  0  dB  or 
greater,  the  advantage  of  the  MLM  scheme  is  obvious  with 
respect  to  the  other  two  functions.  One  also  notes  that  the 
MLM  estimator  is  well  adapted  to  resolving  highly  correlated 
noise  situations;  in  such  cases,  simulations  show  that  the 
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width  of  the  main  lobe  becomes  narrower  but  that  the  mean 
component  of  the  surface  increases. 

B.   MISMATCHING  CASE 

Mismatching  occurs  when  a  parameter  used  in  simulation 
of  the  replica  fields  has  been  incorrectly  estimated  and  is 
different  from  the  one  that  created  the  true  data.  In  the 
eddy  localization  problem,  this  defect  may  be  introduced  in 
several  ways: 

-  a  wrong  measure  of  the  source  frequency, 

-  inaccurate  estimation  of  the  source  or  array  position, 

-  inaccurate  estimation  of  the  source  beamwidth, 

-  insufficient       knowledge       of       the       bottom       loss 
characterization, 

-  oversimplification  of  the  SSP. 

This  list  is  not  exhaustive  and  one  must  keep  in  mind 
that  perfect  matching  almost  never  exists  due  to  the 
impossibility  of  any  acoustic  model  to  solve  the  true 
acoustic  wave  equation  in  the  real  ocean.  Because 
mismatching  prevents  a  close  likelihood  between  the  actual 
and  replica  fields.  It  can  be  thought  of  as  an  additional 
noise  which  the  correlation  matrix  has  not  taken  into 
account.  Mismatching  thus  results  in  a  degraded  estimation 
of  the  total  spatial  correlation  matrix  KT.  The  next 
section  studies  in  greater  detail  two  cases  where 
mismatching  is  created  by  a  change  in  the  bottom  loss 
parameterization  and  where  the  borders  of  the  eddy  are 
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smoother  than  what  we  have  used  previously  to  simulate  the 
real  measured  data. 

1 .   Change  in  Bottom  Loss  Parameterization 

All  simulations  were  conducted  with  the  assumption 
that  the  bottom  was  fully  absorbing.  We  will  now  consider 
the  bottom  to  be  a  perfectly  rigid  surface  with  total 
reflection  and  examine  how  this  new  treatment  modifies  the 
likelihood  functions. 

In  the  absence  of  noise  the  difference  in  phase  at 
each  receiver  of  the  array  for  both  the  perfectly  reflecting 
and  fully  absorbing  bottom  conditions  is  depicted  in  Figure 
4.4.  When  the  bottom  is  treated  as  a  perfect  reflector, 
perfect  matching  will  not  occur  because  all  the  replica 
fields  are  constructed  based  on  the  full  absorption 
assumption. 

Figures  4.5  and  4.6  show  the  result  of  a  simulation 
using  the  MLM  estimator  with  characteristics  r^^  =  10"^^  and 
B  =  1.7  to  represent  the  situations  of  no  mismatching  and 
mismatching,  respectively.  When  mismatching  occurs,  the 
amplitude  of  the  peak  decreases  due  to  an  increase  in  the 
mean  value  of  the  likelihood  function.  The  secondary  lobes 
also  become  larger.  It  is  important  to  note  that  the 
degradation  observed  in  Figure  4.6  does  not  imply  that 
treating  the  bottom  as  a  perfect  reflector  is  less  correct; 
it  only  means  that  the  matching  was  done  improperly;  one 
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must  therefore  be  consistent  in  modeling  the  environmental 
input  variables. 

2 .   Oversimplification  of  the  Actual  Eddy 

The  previously  described  eddy  localization 
simulations  were  run  with  an  excessive  simplification  of  the 
true  medium.  The  simulations  assumed  that  there  were  no 
horizontal  gradients  of  sound  speed  in  and  outside  the  eddy 
boundaries,  i.e.,  the  change  of  SSP  occurred  almost 
instantly  at  the  borders  of  the  eddy.  In  an  attempt  to  be 
more  realistic,  the  next  simulation  was  conducted  after 
adding  four  intermediate  SSPs  between  the  two  previously 
utilized  profiles  (Figure  4.7).  The  first  intermediate  SSP 
was  introduced  4  km  before  the  border  of  the  eddy  and  the 
next  ones  added  every  kilometer  thereafter.  Actual  signal 
values  were  generated  with  this  smoothed  baroclinicity ,  and 
replica  fields  were  computed  as  before  with  the  simplistic 
three-profile  model.  Noise  was  omitted  from  this  simulation 
in  order  to  better  appreciate  the  effect  of  this 
mismatching.  Results  for  the  Bucker  method  are  presented  in 
Figure  4.8.  Comparing  this  plot  with  Figure  3.18  we  see 
that  the  main  peak  decreases  but  localization  of  the  eddy 
remains  possible.  The  implication  of  this  simulation  is 
that  identification  (location)  of  strong  frontal  boundaries, 
such  as  the  north  wall  of  the  Gulf  Stream  or  ice  edge 
fronts,  could  be  fairly  exact  but  weaker,  open  ocean 
mesoscale   eddies  may  pose  more  of  a  difficult  problem 
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(assuming  the  same  number  of  profiles  are  used  to  estimate 
the  replica  fields) . 

The  above  examples  suggest  that  in  cases  where  the 
signal-to-noise  ratio  is  quite  low,  it  is  possible  for 
mismatching  to  hide  the  true  location  of  the  maximum  and  so 
possibly  lead  to  a  failure  of  the  procedure.  As  often  as 
possible,  mismatching  must  be  avoided  by  a  comprehensive 
knowledge  of  the  parameters  used  in  the  replica  fields 
calculations . 

C.   COMMENTS  ON  THE  PROCEDURE 

As  we  were  only  interested  in  using  matched  field 
processing  in  acoustic  tomography,  many  simplifications  have 
been  made  to  run  the  simulations.  Although  the  procedure 
seems  to  be  applicable  and  efficient  in  most  cases,  it  is 
necessary  to  test  its  applicability  to  more  complex 
problems . 

1 .   Environment 

It  was  assumed  in  all  the  simulations  that  only  one 
or  two  parameters  were  unknown,  for  example,  the  surface 
sound  speed  or  the  eddy  location.  For  actual  oceanic 
situations,  many  more  properties  can  be  expected  to  vary  in 
space  and  time.  Extending  this  study  to  a  heterogeneous 
medium  is  theoretically  feasible  but  vastly  increases  the 
number  of  unknowns. 

Modeling  a  shallow  water  environment  has  not  been 
considered   as   a   possible   mechanism   to   speed   up   the 
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calculations  even  though  matched  field  processing  remains 
possible  in  this  kind  of  environment.  Several  limitations 
are  apparent.  A  better  knowledge  of  the  bottom  structure  is 
required.  The  ambient  noise  is  moreover  quite  complex  in 
coastal  waters  and  its  spatial  correlation  matrix  would  be 
difficult  to  model.  As  bottom  interaction  is  important  for 
estimator  calculations  (see  the  mismatching  case) ,  it  is 
possible  to  consider  the  bottom  loss  as  an  unknown  parameter 
and  attempt  to  determine  it  through  matched  field 
processing.  To  examine  this  special  case,  the  replica 
fields  are  generated  using  nine  different  bottom  loss  curves 
and  one  attempts  to  find  the  actual  bottom  loss  curve. 
Figure  4.9  illustrates  the  different  bottom  loss  curves  that 
have  been  used  to  create  the  replica  fields  in  shallow  water 
(300  m)  .  The  maximum  of  the  likelihood  estimator  occurs 
when  the  replica  bottom  corresponds  to  the  actual  loss. 
However,  since  the  various  curves  are  so  similar  in  shape 
and  because  the  bottom  loss  has  only  a  weak  to  moderate 
effect  on  the  transmission  loss,  the  correct  bottom  loss 
curve  is  not  sharply  defined.  This  implies  that  a  correct 
specification  of  bottom  loss  for  a  low  loss  bottom  is  not 
required;  not  so  for  a  high  loss  bottom. 
2 .  Model  Consideration 
a.   Resolution 

In  the  typical  target  location  problem,   the 
medium  parameters  are  generally  considered  constant.   It  is 
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only  necessary  to  run  the  PE  model  once  to  compute  the 
pressure  field  at  different  distances.  However,  when  the 
SSP  becomes  the  unknown  in  the  tomography  problem,  the  PE 
model  must  be  run  for  each  replica.  If  we  wish  to  find  the 
correct  shape  of  the  SSP  from  0  to  300  m  in  the  seasonal 
thermocline  with  a  resolution  of  1  m,  the  model  needs  to  be 
run  10-^^^  times  if  the  sound  speed  at  each  depth  can  have 
ten  different  values.  This  simple  example  illustrates  the 
trade-off  between  resolution  and  computer  time.  In  the 
problem  of  eddy  localization,  where  the  limits  of  the 
perturbation  were  allowed  to  vary  over  11  values,  the  PE 
model  was  run  121  times.  For  a  determination  of  more  than 
two  unknowns,  the  basic  theory  of  the  estimators  is  still 
valid  but  the  representation  of  the  ambiguity  surfaces 
becomes  unachievable  due  to  limitations  in  computer  run 
time . 

b.   Noise  Approximation 

The  correlation  matrix  of  the  ambient  noise  was 
modeled  by  a  symmetric  matrix  decaying  exponentially  around 
the  diagonal.  This  approximation  is  rather  poor,  even  in 
deep  water,  because  the  power  of  the  noise  is  assumed  to  be 
the  same  at  each  hydrophone.  A  better  simulation  would  have 
been  to  consider  the  matrix  of  a  noise  which  is,  like  the 
signal,  a  solution  of  the  wave  equation.  Our  matrix  is 
symmetric  even  though  the  actual  matrix  of  the  complex  noise 
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needs  to  be  hermitian,  because  the  cross-correlations  are 
complex. 

3 .   Correlation  Matrix 

If  the  spatial  correlation  matrix  of  the  noise  were 
known,  it  could  be  introduced  in  the  replica  fields  and  we 
could  deal  with  it  as  with  a  noise-free  problem. 

In  practice,  it  is  not  possible  to  separately 
compute  the  noise  and  the  signal  correlation  matrices.  The 
total  matrix  needs  to  be  estimated  from  the  noisy  signal  at 
the  hydrophones.  Several  techniques  are  available,  as  for 
example  the  Fourier  method  recommended  by  Johnson  (1982). 
As  the  matrix  becomes  only  an  estimate;  we  should  expect  a 
slight  mismatching  and  then  decreasing  efficiency  of  the 
estimators. 
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Figure  4.2   Eddy  Localization.   Bartlett  Estimator. 
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Figure  4.4 
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Hydrophones  for  Fully  Absorbing  ( )  and 
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Figure  4.5   Eddy  Localization.   MLM  Estimator, 
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Intermediate  SSPs  ( )  used  for  Transition 

from  Open  Ocean  SSP  (cl(z))  to  the  Eddy  SSP 
(c2(z)) 
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Figure  4.8   Eddy  Localization.   Bucker  Detection  Factor, 
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tion of  Horizontal  Sound  Speed  Gradients 
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V.   CONCLUSION 

Matched  field  processing  has  been  shown  to  be  an 
efficient  way  to  solve  the  inverse  problem  in  ocean  acoustic 
tomography  when  the  ocean  can  be  characterized  by  a  few 
parameters.  The  estimators  which  have  been  used  (Bucker, 
Bartlett,  Maximum  Likelihood,  etc.)  were  generally  robust 
enough  to  find  the  actual  sound  speed  field  of  the  ocean 
under  usual  conditions  of  noise  power  and  noise  correlation. 

A.   SUMMARY  OF  THE  RESULTS 

For  noise-free  conditions  or  high  signal-to-noise  ratios 
(SNR) ,  all  the  estimators  are  able  to  correctly  determine 
the  actual  unknown  parameter  of  the  medium.  The  Bucker 
detection  factor  was  shown  to  be  the  best  function  when  the 
SNR  was  moderate. 

The  efficiency  of  the  various  estimators,  illustrated  by 
their  spreading  about  the  true  value,  decreased  in  cases  of 
low  SNR  introduced  by  a  large  power  or  a  high  spatial 
correlation  of  the  noise.  The  Maximum  Likelihood  method  was 
shown  to  be  the  best  scheme  when  the  ambient  noise  was 
highly  correlated  because  its  spreading  was  less  sensitive 
to  the  degree  of  spatial  correlation  than  the  other 
estimators. 

The  effect  of  mismatching,  when  introduced  in  the 
simulations,  generated  a  decrease  in  the  efficiency  of  the 
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methods.  Analysis  of  several  degrees  of  mismatching 
indicated  that  unambiguous  results  can  be  expected  from 
matched  field  processing  provided  that  the  parameterization 
of  the  medium  is  exact  enough  to  generate  consistent  replica 
pressure  fields. 

B.   WEAKNESS  OF  THE  SIMULATION 

In  order  to  deal  with  reasonable  computer  times,  only 
the  cases  of  one  or  two  unknown  parameters  were  studied. 
For  the  same  reason,  the  actual  acoustic  pressure  field  was 
compared  with  only  a  few  replica  fields.  This  limitation 
leads  to  moderate  resolution  in  the  results  which  could 
easily  be  improved  by  the  generation  of  more  replica  fields. 

The  weakness  in  modeling  the  noise  field  has  already 
been  emphasized  in  Chapter  III. A.  Further  simulations 
should  be  done  with  a  noise  correlation  matrix,  KN,  that  has 
actually  been  obtained  from  measurements  at  sea.  Further 
work  using  a  more  realistic  parameterization  of  the  ocean 
should  be  done,  e.g.,  with  empirical  orthogonal  functions 
(EOF)  and  using  the  technique  to  estimate  EOF  coefficients. 
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APPENDIX 
FORTRAN  77  PROGRAI^  USED  IN  THE  SIMULATION 


PROGRAM  ESTIMA 
C     ACOUSTIC  TOMOGRAFIIY  USING  MATCHED  FIELD  PROCESSING 
C 

C     THIS  PROGRAM  DRAWS  THREE  KIND?  OP  DETECTION  FACTOR  IN  3D 
C      IT  COMPARES  DUCKER  ,  BAR  TLETP  AND  MI.M  MEIIIODS 
C 

C     THE  NOISE  IS  INTRODUCED  BY  ITS  CORREEAII  f)N  MATRIX 
C 

C     THE  NOISE  MATRIX  IS  PRnpORTlONAl.  10  IDENTITY  MATRIX  IN  CASE  OF 
C      UNCORRELATED  NOISE  BUT  IS  IN  GENERAL  HERMmAN  MATRIX 
C 

C      THE  EXACT(OR  MEASURED)  PHASES  AfJD  MAGNIIUDES  ARE  READ  IN  THE 
C     FILE  EXACT  DATA 

C      THE  COMPUTED  PARAMEIERS  ARE  READ  ON  TliE  FILE  NEAR  DATA 
C 

C      THE  MODEL  ALLOWS  A  MAXIMUM  OF  20  RECEIVERS  DUE  10  THE  USE  OF  THE 
C     PARABOLIC  EQUATION  FROM  THE  EARG  PACKAGE 
C 

C     KS  IS  THE  COMPLEX  MATRIX  OF  MEASURED  PARAMETERS 
C     A  IS  THE  MATRIX  OF  GUESSED  PARAMETERS  WE  WANT  TO  MATCH 
C     NR  IS  THE  NUMBER  OF  RECEIVERS  (tlAXlMUM  20) 

C     NF  IS  THE  COMMON  NUMBER  OF  POSSIBLE  VALUES  FOR  2  UNKNOWNS 
C     THE  NUMBER  OF  FIELDS  WE  tlATCII  IS  ACTUALLY  NF'^NF 
C 
C 

REAL  PPHASE  ( 20 ) ,  ^\\^GU  I  (  20 )  ,  X  ( 1 1 )  ,  Y  (  1  I )  ,  DF  ( 1 1 , 1 1 )  ,  DFCONT (11,11) 

REAL  DART(11,11) .FMLM( 11,11) .MCONFCll . 1 1 ) , BCONT( 1 1 , 11) 

COMPLEX  A (20, 20) .KS(20,20) ,CDF(1J , 1 1 ) , KTf 20 , 20) 

COMPLEX  W(20),CSUtI,r)ETLRfl 

REAL  KN(20,20),NORM 
C     99  DEFINES  THE  FILE  OF  MEASUREDf EXACT)  DATA 
C     98  DEFINES  THE  FILE  OF  DATA  WE  WANT  TO  flATCH 

NR=20 

NF=11  i 

C     READ  THE  MEASURED  VALUES  IN  FILE  EXACT  DATA 

DO  1   1=1, NR 

READ(99,600)  PMAGNI  ( I ) ,  PPHASEd) 

1  CONTINUE 

600   FORMAT(   Ell .4 , 2X,F7 . 4) 

C 

C 

C     COMPUTE  THE  ELEMENTS  OF  MATRIX  KS 

999   DO  2  J=1,NR 

DO  3   K=1,NR 

KS(J,K)  =  PMAGNI(J)^'^PMAGNI(K)-EXP(CMPLX(0.  ,PPHASE(J)- 

PPHASE(K))) 
3        CONTINUE 

2  CONTINUE 
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Vf  iV  Vf  Vf  iV  iV  iV  iV  Vf  Vr  Vf  Vf  iV  iV  iV  iV  iV  ")V  iV  Vr  A  Vr  Vr  iV  iV  Vc  Vr  iV  Vf  iV  Vr  -jV  iV  iV  Vc  Vr  iV  '.V  Vr  Vr  iV  iV  ■jV  Vc  iV  ■>>  iV  Vr  :V  Vf  iV  Vc  Vr  iV  iV  Vr  Vr  Vr  Vr  Vr  ■>V  Vf 

C  DEFINE   THE   CORRELATION   MATRIX   KN   OF  THE   CORRELATED   NOISE 

PRINT'"'        VALUE   OF   SIGflA2    7 
READ^,SIGHA2 

FRINT'^,      VALUE   OF   BETA    7 
READ*, BETA 
DO   A      J-I,NR 
DO   5      1-1, NR 

KN(l,J)=SIGfIA2-EXr(-BETA'''ABS(I-J)) 

5  CONTINUE 
U             CONTINUE 

Vr  Vr  iV  Vr  Vr  Vr  Vr  Vr  Vr  Vr  Vr  Vr  Vr  Vr  Vr  Vr  >V  iV  Vr  Vr  Vr  Vr  Vr  Vr  Vr  Vr  Vr  Vr  Vr  Vr  Vr  Vr  Vr  Vr  Vr  Vr  Vr  Vr  Vr  Vr  Vr  Vr  Vr  Vr  Vr  Vr  Vr  Vr  Vr  Vr  Vr  Vr  Vr  Vr  Vr  Vr  Vr  Vr  Vr  Vr  Vr  Vf 

c 

C  CONFUTE   THE   TOTAL  COVARIANCE   NATRIX    INCLUDING   THE   NOISE 

C  BY   ADDING  THE   NOISE   AND  THE   SIGNAL  CORRELATION   MATRICES 

DO      6   J-l.NR 
DO      7    1=1, NR 

KT(I,J)-KS(I,J)+KN(I,J) 

7  CONTINUE 

6  CONTINUE 
C 

c 

C  BEGINNING   OF   MAIN    LOOP    FOR   IIIE    NF    RUNS 

DO    101      INDEX-] ,NF 
DO    100   JNDEX=1,NF 
C  READ    PARAMETERS    OF    FXPECTPn    FIEI.n    IN   FILE    NEAR    DATA 

DO      8    I-1,NR 

READ(98,600)    PMAGNU  I )  ,  PPHASE  ( I ) 

8  CONTINUE 
C 

C  COMPUTE   ELEMENTS   OF   MATRIX   A 

DO   9      J=1,NR 
DO    10   K-l.NR 

A(J.K}  =  PMAGNI(J)-PtIAGNI(K)-EXP(Cf!rLX(0.  ,    PPHASE  (J) - 

*  PPflASECK))) 

10  CONTINUE 

9  CONTINUE 

CCCCCCCCCCCCCCCCCCCCCCrcCCCCCGC 
C  BUCKER   DETECTION   FACTOR   C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

CDF(INDEX,JNDEX)-(;MrLX(0.  ,0.  ) 
DO    11    J=],NR 
DO    12   K=1,NR 
IF(K.NE.J)    CDF ( INDEX, JNDEX) -CDF (INDEX, JNUEX)+ 

*  A(J,KV''CONJG(KT(J,K)) 

12  CONTINUE 

11  CONTINUE 
C 

C  CDF(--.--)    IS   ACTUALLY   REAL  DUE  TO   PROPERTY   OF   MATRICES 

C  A   AND   KT   HENCE   WE   KEEP   ONLY  THE   REAL   PART  OF    IT 

DF( INDEX, JNDEX )-REAL(CDF( INDEX, JNDEX)) 
C 

C  NORMALIZE  THE   DF  FACTOR 

FACrOR=0. 
FACT=0 . 
DO    13      J=1,NR 
DO    14     K=1,NR 

IF(K.NE.J)FACT0R-FACT0RI-A(J,K)'H:0NJG(A(J,K)) 
I  F(  K  .  NE  .  J)FACT-FACT+KT(  J  ,  K)''C0NJG  (KT(  J  ,  K)  ) 
\k  CONTINUE 

13  CONTINUE 
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FACTOR=SqRT(FACTOR) 
FACT-SQRT(FACT) 

DF(  INDEX,  JNUEX)^OF(  INDEX,  JM[)F,X) /FACT/FACTOR 
nF(INDEX,JNUEX)-AnS(DF(INDEX,jfjl)EX)) 
C 

100  CONTINUE 

101  CONTINUE 
C 

C     NORMALIZE  BY  THE  LARGEST  VALUE 
DMAX=0. 
DO  67  1=1, NF 
DO  68  J=1,NF 

IF  (OF  (I,  J)  .GE.DtlAX)  DNAX-I1F(  I ,  J) 

68  CONTINUE 
67    CONTINUE 

C 

DO  69  1=1, NF 
DO  70  J=1,NF 

DF(l.J)=DF(I,J)/nMAX 
70       CONTINUE 

69  CONTINUE 

C 

C      DISSFLAY  THE  MATRIX   OF  rJORFIALI  ZED  DETECIION  FAC  lOR 

FRINT^'',   DUCKER  FAC  lOR  MA  IRIX 

DO  15   I-1,NF 

WR1TE(6,7  7  7)(DF(I,J),J=1 ,NF) 
77  7      F0RfIAr(lI(2X,F/.  .2)) 

15  CONTINUE 
C 

CCCCCCCCCCCCCCfX'CCCCCCCCfX 
C     DARTLEIT  ESTIMATOR  C 

ccccccccccccuccccccccccccc 

REWIND  98 

C     DEOINNING  OF  MAIN  LOOF  FOR  'l  HE  NF  RUNS 

DO  201   INDEX=] .NF 

DO  200  JNDEX-1 ,NF 
C 
C     READ  PARAMETERS  OF  REPLICA  FIELD 

DO  50  1=1, NR 

READ(96,600j  rMAGNI(I) ,PPHASE(I) 
50    CONTINUE 
C 
C     DETERMINE  THE  NORMALIZED  COMPLEX  VECTOR  W 

NORM=0. 

DO  16  I-1,NR 

N0RM=N0RM  +  PMAGNI(I)^'"'^2 

16  CONTINUE 
NORM=S0RT(NORM) 

C 

DO    17    1=1, NR 

W(l)  =  (i./NORM)*FMAGNI(I)^'^   EXr(CMFLX(0.  ,  PP1IASE(  I )  )  ) 

17  continOe 

C 

C     COMPUTE  THE  DARTLEIT  FACTOR 
CSUM=(0. ,0.) 
DO  18  1=1, NR 
DO  19  J=1,NR 

CSUM=CSUM  +  CONJG  ( W  ( I  )  )  ^'^KT  ( I  ,  J  )  -W  ( J ) 
19       CONTINUE 

18  CONTINUE 

DART( INDEX, JNDEX)=REAL(CSUM) 


c 

200  CONTINUE 

201  CONriNUF. 

c 

C  NORflALIZATION   BY   Till'    HIGHEST  VALUE 

DARHAX-0. 
DO   20      1=1, NF 
DO   21      J=1,NF 

lF(nART(I ,  J)  .GE.DARflAX)    I3ARMAX=DART(1 , J) 

21  CONTINUE 
20    CONTINUE 

C 

DO  22   1=1, NF 
DO  23   J=l ,NF 

DART(I,J)=BART(I,J)/BARMAX 

23  CONTINUE 

22  CONTINUE 
C 

C     DISSPLAY  THE  MATRIX  OF  NORMALI  ZF.D  DETECTION  FACTOR 
^RlNr^   DARTLEIT  FAC'IORS  flAlRIX 
DO  24   1=1, NF 

WRITE (6, 777) (PAR  1(1, J) ,J=1,NF) 

24  CONTINUE 
C 

cccccrccccccccccccccc 

C  HLH   ESTltlATOR   C 

ccccccccccccccccccccc 

REWirn)    9R 

c 

C  INVERT  'IllE   CORRELATION   flATRIX   KT 

CALL  CfMR  I N  ( NR  ,  K  f ,  NR  ,  HE  lERfl ) 
C 
C  BEGINNING   OF   MAIf^    LOOP   FOR   THE   NF   RUNS 

DO    301      INDEX=1 ,NF 

DO  300  JNDEX=1,NF 
C 
C     READ  PARAMETERS  OF  REPLICA  FIELD 

DO  51  1=1, NR 

READ(9B,600)  PMAGNI (I ) , PPHASE(  I ) 
51     CONTINUE 
C 
C     DETERMINE  THE  NORMALIZED  COMPLEX  VECTOR  W 

NORM=0, 

DO  7(-.    I-'  NP 

NC>RM  =  N'ORnf PMAGNI  (I  )^'-->2 
26    CONTINUE 

NORM=S0RT(NORM) 
C 

DO  27  1=1  NR 

W(I)  =  (1./N0RM)^'^PMAGNI(I)-  EXP(CMPLX(0.  ,PPnASE(I))) 
2  7    CONTINUE 
C 

C     COMPUTE  THE  HLM  FACTOR 
CSUM=(0. ,0.) 
DO  28  1=1, NR 
DO  29  J=1,NR 

CSUM=CSUM+CONJG(W(I))>VKT(I,J)-W(J) 
29       CONTINUE 
28    CONTINUE 

FMLM( INDEX, JNDEX)=REAL(1./CSUM) 

300  CONTINUE 

301  CONTINUE 
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c 

C  NORMALIZATION   BY   THE   HIGHEST   VALUE 

Ft]LMAX=0. 
DO   30      1=1, NF 
DO   31      J=1,NF 

IF  (FMLM  ( I ,  J )  .  GE  .  FtlLMAX )    FtlLMAX^FtlLM  ( I ,  J) 

31  CONTINUE 
30    CONTINUE 

C 

DO  32   1=1, NF 
DO  33   J=1,NF 

FMLM ( I , J ) =FM  LM ( I , J ) / FMLMAX 

33  CONTINUE 

32  CONTINUE 
C 

C     DISSFLAY  THE  MATRIX   OF  NORMALIZED  DETECTION  FACTOR 
PRINT*,'  MLM  FACTORS  MATRIX' 
DO  34   1^1, NF 

WRITE(6,7  77)(FMLM(I,J),J=1,NF) 

34  CONTINUE 
C 

c 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C      COMPUTE  THE  SPREADING  FACTOR  OF  THE  ESTIMATORS  C 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
WDF^O. 
WBARr=0. 
WFflLM^n. 
DO  650  1=1 ,NF 
DO  651  J=1,NF 

WDF=WDF  ^  UF(  I ,  n'H'-'^^  • -•  I -'•0  .  )^""-2''' (54  )J-60)''"'^2 
WDAFT=WBART4DARTn  ,J)-(44.4I-40.  )''"-2^>(54t-J-60)-''^2 
WFflLM-WFMLtUFMLfUI  ,Jj'M44.TI-40.  )''^''^2^M54  .  4  J-60  .  )>'" 

651  corrrifJUE 

650   CONTINUE 

c 
c 

PRINT-, '  WDF=  ' ,WDF 
PRINT-. '  WBART=   .WBAPT 
PRINT-,  '  WFMLM=  '  ,WFML[1 


END 
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